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Abstract 


The Semi-Discrete Galerkin Finite Element Modelling of 
Compressible Viscous Flow Past an Airfoil 

by 

Andrew J. Meade, Jr. 


A method is developed to solve the two-dimensional, steady, compressible, turbulent 
boundary-layer equations and is coupled to an existing Euler solver for attached tran- 
sonic airfoil analysis problems. The boundary-layer formulation utilizes the semi-discrete 
Galerkin (SDG) method to model the spatial variable normal to the surface with linear finite 
elements and the time-like variable with finite differences. A Dorodnitsyn transformed 
system of equations is used to bound the infinite spatial domain thereby permitting the use 
of a uniform finite element grid which provides high resolution near the wall and automati- 
cally follows boundary-layer growth. The second-order accurate Crank-Nicholson scheme 
is applied along with a linearization method to take advantage of the parabolic nature of the 
boundary-layer equations and generate a non-iterative marching routine. The SDG code 
can be applied to any smoothly-connected airfoil shape without modification and can be 
coupled to any inviscid flow solver. In this analysis, a direct viscous-inviscid interaction 
is accomplished between the Euler and boundary-layer codes through the application of a 
transpiration velocity boundary condition. Results are presented for compressible turbu- 
lent flow past NACA 0012 and RAE 2822 airfoils at various freestream Mach numbers, 
Reynolds numbers, and angles of attack. All results show good agreement with experi- 
ment, and the coupled code has proven to be a computationally-efficient and accurate airfoil 

analysis tool. 
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Chapter 1 
Introduction 


The successful application of computational fluid dynamics to the design and analysis 
of two-dimensional airfoils in transonic flow regimes has been accomplished by a large 
number of researchers in the past fifteen years. Analysis of the entire flowfield around 
wing sections can generally be performed through the use of two main techniques. The 
first technique is to solve some practical form of the Navier-Stokes equations for the entire 
flowfield. The second technique involves the solution of the boundary-layer equations 
in the viscous flowfield region and either the potential or Euler equations in the mviscid 
flowfield region. Iteratively solving the viscous and inviscid equations while enforcing a 
compatibility condition then yields a solution for the entire flowfield. 

The Navier-Stokes methods, while modelling most of the physical flow mechanisms 
and providing accurate results, require large amounts of computer time and storage. The 
coupled viscous-inviscid methods are generally 30-500 times faster than the Navier-Stokes 
methods and generate results of adequate accuracy [33]. Advances in computer hardware 
technology are constantly narrowing the computational advantage held by coupled methods 
over the Navier-Stokes methods. However, coupled methods should continue to be partic- 
ularly important in an interactive airfoil design environment where near real-time flowfield 
analysis is desired. 

The numerical solution of the classical boundary-layer equations has traditionally been 
accomplished through the use of finite differences [1]. The finite element method has been 
used only since 1972 to obtain numerical boundary-layer solutions, even though the method 
itself has been in existence since 1915 [11]. Traditionally, the finite element treatment of 
two-dimensional boundary-layer flow has involved the use of finite differences in the 
stream wise direction to capitalize on the parabolic nature of the boundary-layer equations. 
Coordinate transformations, which reduce the number of dependent variables or generate 
computational grids suited for finite elements, have also been used. 

In 1960, A.A. Dorodnitsyn developed and applied a set of boundary-layer transfor- 
mations which are well-adapted to the finite element method of solution [15]. Initially, 
the transformations were used along with the method of integral relations to solve various 
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classes of supersonic boundary-layer flows [2]. Fletcher and Fleet successfully applied the 
finite element method to the laminar and turbulent incompressible Dorodnitsyn form of the 
boundary-layer equations [3], [18]. Meade and Strong have extended the method to solve 
laminar compressible flow about cones and airfoils, respectively [17], [4]. 

Prandtl proposed that the inviscid pressure distribution could be determined to a higher- 
order accuracy by recalculating the potential flow while accounting for the displacement 
thickness of the boundary layer [5]. Perhaps the most notable transonic airfoil analysis 
code, employing Prandtl’s direct interaction procedure to couple Green’s [6] lag-entrainment 
boundary-layer code with an inviscid full-potential code, is the viscous Garabedian and Korn 
program developed by Collyer and Lock [7]. Another notable example is the GRUMFOIL 
program developed by Melnik, Chow, Mead, and Jameson [8]. A great deal of research has 
been performed to develop viscous-inviscid coupling mechanisms which do not have the 
disadvantages associated with the direct displacement thickness approach. The transpiration 
velocity boundary condition, as suggested by Lighthill [24], has been used as the viscous- 
inviscid coupling mechanism by Van Dalsem, Steger, and Rao with success [26]. 

The viscous, compressible, transonic airfoil analysis method presented in this text is 
based on the direct viscous-inviscid interaction of finite element boundary layer and fi- 
nite difference Euler codes. The present work extends the semi-discrete Galerkin (SDG) 
boundary-layer formulation to include turbulent flow. A direct transpiration velocity 
viscous-inviscid interaction approach will be used to couple the SDG method with an 
innovative Euler solver (GAUSS2) which employs a shock fitting technique. The coupled 
codes will be used to analyze flow about a NACA 0012 and RAE 2822 airfoil for attached, 
turbulent, compressible flow. 
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Chapter 2 

Finite Element Method 


The finite element method is a numerical analysis technique for obtaining piecewise approx- 
imate solutions to the governing equations of a wide variety of engineering problems. The 
principle of the finite element method is to replace a continuum, having an infinite number 
of unknowns, with a discretized domain of assembled elements, having a finite number of 
unknowns. The unknown field variable is expressed in terms of assumed approximating 
functions within each element. The approximating, or interpolation, functions are defined 
in terms of field variable values at specific points or nodes. Each element has a prescribed 
number of nodes which may be on the boundary, where connections to other elements are 
made, or in the interior of the element. Thus, the nodal values of the field variable and 
the interpolation functions completely define the behavior of the field variable within the 
elements. 

The solution to any continuum problem by the finite element method is accomplished 
in the following steps [10]: 

• Discretize the continuum 

• Select interpolation functions 

• Find the element properties 

• Assemble the element properties to obtain the system of equations 

• Solve the system equations 

2.1 Finite Element Approximation 

Consider <j>(x) as an approximate, or trial, solution to the one-dimensional field variable 
<f>(x), which can be written as 

M 

<f>{ x ) = 

1=1 


( 2 . 1 ) 



4 


where 4>i are the nodal unknowns, N,(x) are the interpolation functions, and M is the total 
number of nodes or nodal unknowns. The derivative of 4>(x) is approximated in this finite 
element representation by 

d$“dNi 
dx Til dx 


■4>i- 


( 2 . 2 ) 


Referring to Figure 2.1, it can be seen that when linear interpolation functions Ni are used, 
4>{x) interpolates the function <f>(x) linearly over each element. 



Figure 2.1 Finite Element Representation of <f>(x) 


2.2 Interpolation Functions 

Interpolation functions are normally chosen to be locally defined polynomials within each 
element. It can be seen from Figure 2.2 that linear interpolation functions fall from a 
maximum value of one at a particular node to zero at the two neighboring nodes and 
are zero throughout the rest of the domain. Therefore, even though Equation 2.1 is a 
global equation, only two nodal unknowns and two interpolation functions make a nonzero 
contribution in any particular element. As shown in Figure 2.2, the local shape functions 
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satisfy the following conditions in each element (e): 

jvf e) (x) = 0 if x not in element (e) 

M 

Y N- e \x) = 1 for all x in element (e) 

.=i 

It is also necessary that the interpolation functions be chosen in such a way that the field 



Figure 2.2 


Linear Interpolation Function N- e \x) 


variable <f>(x) and any of its derivatives, up to one order less than the highest order derivative 
appearing in the weak form of the equation, be continuous at the element boundaries [11]. 
The linear interpolation functions can be determined by using Lagrange polynomials and, 
for elements A and B, take the form 


N, b 


X X i 

ar,_i - Xi 

X 3't + l 


N t A = 


N , B + . = 


X — £t-l 
X — X{ 


^t+1 Xi 


X{ ^ t -h 1 


(2.3) 
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2.3 Method of Weighted Residuals 

The method of weighted residuals is a technique for obtaining approximate solutions to 
partial differential equations. The method of weighted residuals is one of many approaches, 
namely the direct approach, the variational approach, or the energy-balance approach, used 
to determine the finite element matrix equations which express the properties of individual 

elements. 

Applying the method of weighted residuals involves basically two steps. The first 
step is to assume the general functional behavior of the dependent field variable so that the 
given differential equation and boundary conditions are approximately satisfied. A residual, 
which is required to vanish over the entire solution domain, results when the approximation 
is substituted into the original differential equation and boundary conditions. The second 
step is to solve the equation(s) resulting from step one for a particular function. 

Consider finding an approximate functional representation for the field variable <j> which 
is governed by the differential equation 

L(<t>) - f = 0, ( 2 - 4 ) 

where L represents the differential operator and / is a known function of the independent 
variables. The differential equation resides in the domain D bounded by the surface S, 
where proper boundary conditions are prescribed. The unknown exact solution for 4> can 

V 

be approximated by 4> as 

Af 

<f> ~ ^ Nifa, 

1 = 1 

where JV, are the assumed functions and <f>i are the unknown parameters. The M functions 
Ni are usually chosen to satisfy the global boundary conditions. 

When 4> is substituted into Equation 2.4, a residual R results from the approximation 

and is given by 

R = L{$) - f. ( 2 - 6 ) 

The method of weighted residuals seeks to determine the M unknowns 4>i in such a way that 
the error R over the entire solution domain is minimized. This is accomplished by forming 
a weighted average of the error and specifying that this weighted average vanish over the 
solution domain. Therefore, choosing M linearly independent weighting functions Wj and 

insisting that if 

J W j [L(i)-f}dD = WjRdD = 0, j = 1,2, ..., M 


(2.7) 
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then fi«0. The form of error distribution principle used in Equation 2.7 is dependent 
on the choice of weighting function. There are a variety of weighted-residual techniques 
which can be employed; the most popular error distribution principle is the Galerkin 
Method. The Bubnov-Galerkin (classical Galerkin) method uses the interpolation functions 
Nj as the weighting functions Wj, while the Petrov-Galerkin method specifies W 3 as some 
modification of N r 

2.3.1 Petrov-Galerkin Method 

The Petrov-Galerkin method is used in applications of the method of weighted residuals 
when specific requirements must be imposed on the finite element solution. The weighting 
function in this method is represented by 

= Pj, ( 2 - 8 ) 

where Pj is an analytic function similar to the Bubnov-Galerkin interpolation function 
but with additional terms or factors to impose the specific solution requirements. The use 
of the Petrov-Galerkin method has been based in part on its ability to produce asymmetric 
weighting functions which force diagonal dominance of the finite element matrix equations 
and reduce the oscillatory solution behavior in convection-dominated fluid flows [11]. 


2.4 Semi-Discrete Galerkin Method 

The semi-discrete Galerkin method is a hybrid finite element and finite difference numerical- 
analysis technique which uses a finite element representation for the spatial variables and 
models the time or time-like variables by finite differences. The semi-discrete Galerkin 
(SDG) method is best demonstrated by its application to the one-dimensional, unsteady, 
nondimensionalized Heat Conduction equation: 

5^ _ = o (2.9) 

dt dx 2 

The initial and boundary conditions for <f>(x, t), over the interval 0 < x < 1, are given by 

<f>(x,0) = <f> 0 (x), 4>{0,t) = 0 , <HM) = 1- (2-10) 

Substituting an approximate solution for <f>, 

M 

fa,*) = H N.-(z)<M0> 

t=l 


( 2 . 11 ) 
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into Equation 2.9 and applying the method of weighted residuals gives 

f 1 n 


[ w ‘ (I - 


dx = 0, 


( 2 . 12 ) 


where j = 1 , 2, M. The classical Galerkin method can be used by taking the weighting 
function as the interpolation function: 


Wj = 


Therefore, Equation 2.12 becomes 




(2.13) 


(2.14) 


In order to satisfy the continuity requirements for linear interpolation functions, it is nec- 
essary to reduce the second partial differential by applying the Green-Gauss theorem. The 
Green-Gauss theorem in one-dimension is simply an integration by parts which gives 


A dl , f 1 dNj d(j> , . r d]> 

L N 4 iz --L^si dx+N ‘ai 


(2.15) 


The last term in Equation 2.15 represents the natural boundary conditions. If Neumann 
boundary conditions are present, it would be appropriate to replace the boundary term with 
the given condition. However, Dirichlet boundary conditions are specified so the bound- 
ary term remains incorporated in the governing equation. Substituting the approximate 
solutions for 6 and §| into Equation 2.15 produces the following finite element equation: 


A ^ dh , f 1 dNj ^dN iA , , ^ dN, 

l N E«^-4x = -/ o — g + Nj E dx * o 

Rewriting in a more convenient form by introducing the coefficients, 

t\ A dNj dNi 

= l NjN'dx, A2„ = l 


( 2 . 16 ) 


gives 


M ji M M /IN 

X>',f —£** + *£%* • < 217 > 

i=l ut i=l «=i o 

A time-independent discretization in x (uniform) allows the solution of Al and >42 to be 
obtained prior to the rest of the the solution process. The elements of matrices Al and 
A2 can be solved exactly by Gaussian quadrature if lower-order polynomials are chosen as 
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interpolation functions. Also, the matrices become tridiagonal if the interpolation functions 
are chosen to be linear. Modifying A2 to absorb the natural boundary conditions and 
renaming as A3 result in the following system of ODE s: 


,= 1 al i = 1 


(2.18) 


The boundary conditions in nodal form are given by 


<f> 1 =0, <f>M — 1- 


To utilize the parabolic nature of the governing equations, a marching routine is invoked 
by using a finite difference discretization for 


4>? +l - 

At 


(2.19) 


t n +l _ t n 

where n denotes the time level. The resulting A# 1 * 1 may then be solved for using the theta 
method: M r M m 

Y'AIjM ?* 1 =-A< 


( 2 . 20 ) 


0'£A3 JI <f>? + '+(l-0)E A3 ^ 

. t=l 

The value of 0 controls the degree of implicitness. A number of schemes which depend on 
the value of theta may be employed: 


0 = 0.0 fully-explicit Euler forward 

0=1.0 fully-implicit Euler backward 

0 = 0.5 second-order accurate Crank-Nicholson 

Choosing the Crank-Nicholson technique and linearizing Equation 2.20 converts the system 
of ordinary differential equations into a system of linear algebraic equations given by 

A/ M 

Y,(A\ ji + 0AtA3 ji )Atf +1 =-Af£A3 Jt #\ (2.21) 

,=i 

A detailed discussion of the linearization is given in Section 3.3. 

The right-hand side of Equation 2.21 is known for a given time increment A t and initial 
condition #*. The left-hand side array forms a tridiagonal matrix of size ( M x M) that 
may be efficiently solved by the Thomas Algorithm which is ideally suited for equations of 
this type [13]. The Thomas Algorithm requires an order of M operations ( 0(M )) which 
is more efficient than the 0(M 2 ) operations necessary for Gaussian elimination. After the 
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system of equations has been solved, A#* +1 is added to the known value to give #* +1 . 
The marching routine continues as ^” +1 is used in the right-hand side of the equation to 
solve for A <£" +1 at the next time step. The time step is varied to obtain a desired accuracy, 
instead of performing iterations, which results in a computationally efficient algorithm. 


2.5 Group Finite Element Method 

Finite element treatment of the nonlinear convective terms, which are present in most 
flow problems, is traditionally accomplished by the introduction of a separate approximate 
solution for each contributing variable in the nonlinear terms. The group finite element 
method permits the nonlinear convective terms to be represented without introducing a 
separate approximate solution for each variable. Since the added connectivity of separate 
approximate solutions is avoided, the group finite element formulation can model nonlinear 
terms and avoid having products of nodal values over all connected nodes in a particular 
element [12]. 

The group finite element method consists of transforming any convective terms into 
a divergence form and then employing supplementary solutions for these terms. The 
group method is best demonstrated by its application to the one-dimensional, unsteady, 
nondimensionalized Burger’s equation: 


d<j> d<j> 1 &<f> 
~dt +<t> dx~ Re dx 2 

Transforming the convective term into a divergence form, 

,<H _ 1 d(4> 2 ) 

* dx ~ 2 dx ' 


( 2 . 22 ) 


(2.23) 


and substituting into Equation 2.22 gives 

d<f> \ d{£) 1 &4> 

dt + 2 dx Re dx 2 

The supplemental approximate solution is then introduced for the group <f> 2 as 


(2.24) 


M 

$ 2 (M) = ^L, Ni(x)<f> 2 i(t). 

t=l 


(2.25) 


After substituting the supplemental approximate solution into Equation 2.24 and applying 
Galerkin’s method with linear elements on a uniform grid, the following system of ODE’s 
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is produced: 


d_ 

dt 



2 1 \ 1 , , , . <£.+i - <£.- 1 

+ 3 ^' + 6^ ,+1 ) + 2 + ^ <+1 2Ax 

1 / <£■- i ~ 2<£, + <£i+i \ 

_ Re V Ax 2 / 


(2.26) 


The traditional finite element formulation produced the following similar system of ODE s. 


d_ 

dt 



2 1 . 

+ 6^ ,+1 


^ + X (<£i-l + <£• + &+l) 

= w*zl 

Re V 


<£ i + 1 _ ^*-1 


2Ax 

_i - 2<£, + <£.+i 


) 


(2.27) 


As seen in the second term of Equations 2.26 and 2.27, the group method produces 
a computationally more economical finite element form of Burger’s equation by reducing 
the nodal connectivity of the convective term. The group finite element method becomes 
progressively more economical as the order of nonlinearity or the number of dimensions 
increases and generally produces a more accurate finite element scheme [12]. 
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Chapter 3 

Boundary Layer Formulation 

The boundary-layer approximation of the Navier-Stokes equations is valid if the viscous 
flow region, prior to separation from a body, is assumed to be thin. That is, the boundary- 
layer thickness is much smaller than the characteristic length of the body in question. 
Prandtl proposed that the following assumptions could be made about thin shear layer flows 

[14]: 

• negligible body forces 

• negligible normal viscous stresses 

• negligible normal pressure gradient 

• normal velocity <C tangential velocity 

• tangential velocity gradients <C normal velocity gradients 


3.1 Governing Equations 


The equations 
given in terms 
Continuity: 


of motion for steady, compressible, and turbulent boundary-layer flow are 
of dimensional variables: 


m (?5) + Ty W = 0 


(3.1) 


Momentum: 


?s a? + p % 


dp <>_ 

dx dy 


' _ . du 

+ 


(3.2) 


Energy: 

„dR ~dH 

pn Si + pv ls 


j_a_ 

Pr dy 


an l /, M9L ± .,i (±\ 

^ + p,) ~di j + V 1 " Pr) as r as U J. 


(3.3) 
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The Prandtl number is assumed constant, and the perfect gas assumption serves as the 
equation of state of the fluid. The gross effects of turbulence in the boundary layer are 
accounted for using the eddy- viscosity model 


- du 

PtTTZ = -pu'v’, 
ay 


(3.4) 


where the tilde denotes average values and the prime denotes fluctuating values [9]. 
Dropping the tilde and prime notation, the dimensional velocity components are denoted 
as u in the local x coordinate and v in the local y coordinate. The dimensional pressure, 
density, laminar viscosity, turbulent viscosity, and total enthalpy are shown as p , p , p, pt> 
and H , respectively. 

The initial conditions for velocity and temperature are provided at x = x 0 and repre- 
sented by 


u(x 0 ,y) = u 0 {y), v{x 0 ,y) = v 0 {y), and T(x 0 ,y) = T 0 (y), 


(3.5) 


where u 0 (y), and T 0 (y) are known quantities. 

Boundary conditions are prescribed at y = 0 and y = oo. The boundary conditions 
for velocity are obtained by noting the no-slip condition at the surface of the body and that 
the tangential velocity approaches the magnitude of the inviscid velocity at the edge of the 
boundary layer. The velocity boundary conditions are given by 

u(x,0) = u(f ,0) = 0 and u(x, cx>) = u e (x), (3.6) 

where u e {x ) is the velocity at the edge of the boundary layer obtained from a solution of 
the inviscid equations of motion. The heat-transfer boundary conditions for temperature 
are given by 

f(x,0) = T w (x) and T (x, oo) = f e (x), (3.7) 

where w and e denote wall and boundary-layer edge conditions, respectively. If an adiabatic 
wall is prescribed, the following boundary conditions for temperature are required: 

= 0 and T(x, oo) = T e (x) (3.8) 

v=o 

The equations of motion. Equations 3. 1-3. 3, can be nondimensionalized by employing 
a characteristic length L, critical speed of sound a m , and free-stream density poo • The 


dT 

dx 
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dimensionless variables are defined as follows: 


x = 


V = 


u 

a* 


v 

v = — 
a* 


P — T 


P P 

P fi _ . T 

Or^ Poet a L* 


„ +l(fi 2 + 5 2 ) 
H= a - 2 


(3.9) 


Poo Pood 

After substitution of the nondimensionalized parameters and simplifying, the continuity 


equation becomes 

Jr (pu) + - (pv) = 0. (3.10) 

ax at/ 

Since pressure is only a function of x and is transmitted unchanged through the boundary 
layer. 


dp _ dp _ dp e 
dx dx dx ’ 

where p e is the dimensionless pressure at the edge of the boundary layer. Substituting the 
pressure term along with the nondimensionalized parameters into the momentum equation 


gives 


du 


du 


dpe , 9 


Tx + fV Ty = ~~dx + dy 


du 

O' + V.) 


(3.11) 


Similarly, the energy equation takes the form 


dH dH 

pu -^ +im i; 


\_d_ 

Pr dy 


v 9H' 

(P + Pt)~^ 



(p + Pt) 


d_ 

dy 



(3.12) 


3.2 Traditional Dorodnitsyn Formulation 


A. A. Dorodnitsyn applied the method of integral relations to the compressible boundary- 
layer equations and transformed them into a form resembling the incompressible equations 
[15]. Dorodnitsyn used the following expressions to transform the body-normal x and y 
coordinate system to a £ and rj computational plane and smooth both velocity and density 


over the boundary layer: 


The Dorodnitsyn formulation 



1 

f u e p e dx 

^ooPoo J 

'o 

1) = 

U e 

7 Jo pdy 


(UooPoo 

) 2 Jo 


utilized a normalized velocity u 


(3.13) 

(3.14) 


U 



(3.15) 
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and an intermediate variable w given by 

„ * UooPoo dt) 


. “OO^OO ''• 1,1 1 

W = U -T 1- I rp ■ 

U e p e OX \Po o) u c * 


(3.16) 


Employing the Dorodnitsyn relations results in the following form of the continuity equa- 


du dw _ 
d( + dr) 


(3.17) 


It is assumed that the variation of laminar viscosity across the boundary layer can be 
represented by the linear relation 

-S- = , <3- 18 > 

Li d oo 


where C 0 is the Chapman-Rubesin constant. A full explanation of the assumption is given 
in Reference [16]. Applying the Chapman-Rubesin relation, ideal-gas law, Dorodnitsyn 
relations, and a viscosity parameter a 

<7 = 1 + — , ( 3 - 19 ) 


results in the following form of the momentum equation: 


6 ?H + *2“ = I^(l -an+CoAf^) (3.20) 

di dr) u e d£ v ' dr)\ dr)) 

Using a new variable 5 which relates enthalpy and the ratio of specific heats (7 = c p /c v ) 


- > - (4) ■ 

with the Dorodnitsyn relations results in the following form of the energy equation. 


~ds ds 

U 3( +U V 2 


= «A+2g, 

u e d£ Prdr) \ dt) ) 


(7 + 1) \Pr ) dr) \ dr) 


(3.21) 


(3.22) 


The Dorodnitsyn formulation, representing the two-dimensional, compressible, boundary- 
layer equations in incompressible form with respect to the independent variables, t and 7, 
is therefore given by Equations 3.17, 3.20, and 3.22. An integral form of these equations 
can be obtained by applying the method of weighted residuals. 

After using a general weighting function f{u) and summing the products of Equation 
3.17 x / and Equation 3.20 x % the first Dorodnitsyn equation appears as 


«M + ^> = ±±i(i-a *) £ + \ a ir 

dt + dr) Ue dt ' ' du du dr) \ dr) 


(3.23) 
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Integrating the above equation with respect to 77 and reducing gives the following: 

± f (3 ' 24) 

Applying the Crocco transformation to the integral equation changes the independent vari- 
able 77 to u [17]. The transformation makes use of the nondimensionalized shear stress and 

isgivenby r nw 

u — I rdr], (3.25) 

Jo 


or 


r = 


du 

di)' 


(3.26) 


Substituting the Crocco transformation into Equation 3.24 and noting that the limits of 
integration are changed from tj = 0— » 00 to u = 0 — ♦ 1, gives 

i + */ii - £ (> *) + c °£ ss <"> "• (327) 

Assuming that there is no surface injection, v, and hence w, is zero at the wall. Therefore, 
if f(u) is chosen to vanish at the edge of the boundary layer, the traditional Dorodnitsyn 
formulation for the first integral equation reduces to 

— [' uf -du = — — /' (1 -u 2 )^z-du + Co f ^-^((TT)du. (3.28) 

d£ Jo t u e d£ Jo ' ' du t Jo du du 

After using a general weighting function /(«) and summing the products of Equation 

3.17 x 5 x /, Equation 3.20 x 5 x and Equation 3.22 x /, the second Dorodnitsyn 

equation appears as 

df 


d ( ufs ) d ( wfs ) 


di 


dr] 


_ 1 du e „ 1 du e / 2 \ 

-2 -u/s + — — ) 


li e 


li e di 


du 


+ 


+ 2 Co 




(7 - 1) 


(7 + 1) 

Integrating the above equation with respect to 77 and reducing gives the following. 


(3.29) 


-jr [ ufsdi] + wfs\™ 
df Jo 


= -2 


1 du e 


f°° 

/ iifsdrj 
Jo 


+ 


1 du e 


+ 


/ (l - u 2 ) — sify 

li e di JO ' U e di Jo V ' 

df d ( du\ . Co [°° f d ( ds 
Co Jo dZdv ^ + Pr h dr > 


di) 


dr] 
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Using the Green-Gauss theorem to reduce the second-order derivative of s with respect to 
, h and again applying the Crocco transformation while noting the previous assumptions for 
/(u) and w, gives the modified Dorodnitsyn formulation for the second integral equation: 


-t 

d£ Jo 


f 1 S , _ 1 du e f l „ f s 1 du e y* / _ f2 \ £ JA 

l Vr iu = ~ 2 7 ,liJo uf T du+ U' a Jo >dur 


r \ df d Co r 9s 1 Co [ l df ds 

+ Co /o Tu S Tu ( ° T)iiJr Tr ! ° T Tu-TrL du” d 6 


+ 2C "77Tt1 4 - 0 Jo ^ 


(3.31) 


T ‘ u ( 7 + l) VPr / 7o ' 

In order to permit a more physically-oriented solution procedure, it is desirable to use 

the independent variable x in place of £. Making use of the fact that 

d d dx , dx UooPoo 

— = 7 — . where -77 = > 

d£ dx dC d£ u e p e 

Equations 3.28 and 3.31 become 

± t'uf = f (l - fi ! ) %>-d& + Co-^- [' %r-^(°T)du, (3.32) 

dx Jo T tie dx Jo ' > dux UooPoo^O du 


-I' 

dx Jo 


I'uf-du = -2- d -^ !' 

Jo T U e dx Jo T 


. [' h-&)£idu + a-^ [' ££(.<rr)d& 

U e dx Jo ' 2 t(tl T U00P00 do du du 


Co u e p e 5s Co UePe f !!L aT !li l 
+ PrUoopJ 07 du 0 Pru^Pook du du 


+ 2Co 
+ 2Co 


(7 + 1) 
(7 ~ *) 
(7 + 1) 


4-0 

! u «Pe f 

UqoPoo Jo 

( 1 _l) 

, «eP e r 

\Pr 1 J 

^oo Poo •'0 


(3.33) 


respectively. 

The boundary conditions for t and s in terms of the Dorodnitsyn formulation are given 
by 

t( 1) = 0, s(0) = s. = (l - jr) and a(l) = 0, (3.34) 
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in the case of heat transfer at the wall. For an adiabatic wall, the boundary conditions for 6 
are specified as 

= 0 and 5(1) = 0. (3.35) 

u=0 

Numerical experiments have shown that the solution for r was generally inferior when a 
wall boundary condition was imposed; therefore, the wall boundary condition available 
from Equation 3.2 is not used [18]. 

It has been shown that the Dorodnitsyn formulation reduces the nonlinear partial dif- 
ferential equations governing two-dimensional, compressible, and turbulent flow to a set 
of uncoupled integral equations by seeking the proper weighted combinations of variables 
and equations. The integral equations are explicitly independent of density and represented 
in the independent variables x and u. The transformation from a body-normal coordinate 
system ( x and y) to a computational plane (x and u ) is depicted in Figures 3.1 and 3.2. One 
benefit of this coordinate transformation is that the infinite domain in the y direction has 
been replaced by a finite domain in the u direction. A greater benefit is that the uniform 
grid in u automatically captures downstream boundary-layer growth and spatially provides 
high resolution near the wall. 


ds 

du 



3.3 Semi-Discrete Galerkin Formulation 

In order to apply the semi-discrete Galerkin method to the integral Dorodnitsyn equations 
developed above, the general weighting function f(u) must be replaced by a set of linearly 
independent functions fj{u). The weighting function is introduced as 

/,•(«) = (1 - fi)ty(fi), 3 = 1 , 2 ,..., A/ 


(3.36) 
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Figure 3.2 Computational Plane 


where N } (u) is the linear interpolation function at a particular node j, and M is the total 
number of nodes. The Petrov-Galerkin finite element method is utilized by introducing 
the term (1 - u) to satisfy the requirement that fj(u) equals zero at the outer edge of the 
boundary layer. The trial solutions for the dependent variables r, s, a , and or are 
introduced using the group finite element method and take the following form. 


1 

T 


1 

(1 - u) 


M 


t=i 


(3.37) 



(3.38) 




20 


s = (l-u)'£ l N i (ti)b3 i (x) 
«= 1 

s M 

- =^Ni(u)b4,{x) 

T U 

M 

a = Ni{u)b5i{x) 

1=1 

M 

< 7 T = (l-fi)]TA , ,(8)ft6,(x) 


(3.40) 


(3.41) 


(3.42) 


1—1 

Substituting the approximate solutions and the weighting function into Equation 3.32 
produces the following finite element equation: 

, i i M 


+ Co 


HqcPoo ^ 


/;[d- «)§ - «*] [(* - «>g f < - 1 H <3 - 43) 


Removing all terms that are independent of u from the integral simplifies the Equation 3.43 

as follows: 

First term: 


±r 

dx Jo 


1 1 M 1 w ,1 

l fi[(i = iN,Nii " 


Second term: 




Third term: 


^ooPoo •*'0 




M dN 

v- /' f„ 

Si, (1 - U) 1T 


= Co^Ei (i-«)^-*<«“* 

UooPoo "1^0 L du J L ° U J 
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Therefore, after introducing 


I<\ji = j uNjNidu, 
A'3j, = /'(l +«)[(!- ^ 


du, 


Ci = 


1 du e 


and 


C 2 = Co 


ti e dx ’ 

UePe 


(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 


UooPo 


into Equation 3.43, the following system of first-order ordinary differential equations is 
produced: 

A/ //i A/ A/ 

V A'lji Kiji Hi + C 2 £ K2„ b6, (3.49) 

i-1 dx ,=i »=i 

After substituting the approximate solutions into Equation 3.33, removing all terms that 
are independent of u from the integral, and introducing 

KAji = [ u(l - u)NjNidu, 


I<5,, = J\l - & 2 ) (1-u) 




du 


N{du , 


I<6ji = f (1 - ufNjNidu , 

J 0 

r \ r, „,dAr, »/ 

K 7ii = / o fi(i-«w[d-ii)^--JV. 

*«** = /’(>-«) [(> - - N ] [(l-fi)^-^]^’ 

K9 jik = f\ 1-fi) 

0 






du 


A'lOj.jt = (1 - fi) 2 ^ 


ft = ft K , 


du 

L 

(1 -t l)^-N k 
du 

1 


(3.50) 

(3.51) 

(3.52) 

(3.53) 

(3.54) 

(3.55) 

(3.56) 


(3.57) 
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and 


a -«*{$!(£ -O' 


(3.58) 


the following system of first-order ordinary differential equations is produced: 


M 


M 


M 


£ A’4ji ^ = -2C.EA-4,, M, + C,X>'5„ M, 

1=1 


t=l 

M M 


t=l 


M M 


+ C 2 ry] AT 8 j,ib b6ib3 k + Ci 53 13 ^ 10 i‘* 

,=1 Jt=l «=1 k=\ 

MM W 

— C 3 ^ K9jik 66,63^ + C 4 ^ 


t =i i=i 

M 

+ C 4 ^/n ii 66, 

t=l 


(3.59) 


The Galerkin integrals (A r 1 , ..., A" 10) are further defined in Appendix A. 

The second-order Crank-Nicholson scheme, obtained by using 9 = 0.5, produces an 
efficient implicit algorithm for marching the solution in the x direction. Equation 3.49 is 
approximated by 


m r i 

£ I<\ jt A61" +1 = Ax [OF If 1 + (1 - 9)Fl ]\ , 

.=1 


where 


and 


M *’• 

F\j = Ci ^ A'3ji feli + C 2 A2jt fe2,, 

*=i 


M 

£ 

i = l 


A61" +1 = 61” +1 - 61?. 


(3.60) 

(3.61) 

(3.62) 


The superscript n denotes a particular streamwise, or x, position. 

In order to develop a non-iterative marching routine, it is necessary to linearize AT" + 
by an expansion about the nth level following the approach of Briley and MacDonald [19]: 


fi 7 ‘ = ri > mS:) 


(3.63) 


Substitution of the linearization into Equation 3.60 produces the following system of equa- 
tions in 61: 

YKTjitib 1" +1 = FT, 


(3.64) 



23 


The variables I\T and FT are defined as 

M M 


i= 1 

- c, n+1 


52 KTji = E K1 H - eAx [ C " +1 2 


and 


{ Afl 

0C, n+1 + (l - ^c?] E^ 3 * 61 " 

A/ 

Similarly, applying the Crank-Nicholson scheme to Equation 3.59 gives 

M . T 

55 K4ji A64" +1 = Ax [6F2] +l + (1 - ^)F2"j , 


where 


M 


M 

F2j = -2C\52 I<4jitei + c \52 K5 i iM ' 

i = l *=1 

M M 


M M 


+ C 2 E E K%jik 66,63* + Ci 5252 KlO jik 66, A3* 

1=1 fc=l < =l fc =l 

M M M 

— C$ E 55 C9jik 66,63* + Ca, 55 A 


l=\ k = 1 

M 

+ Ct'EMiiMi , 

t=l 


1=1 


and 


A64” +1 = M” +1 - 64". 

Performing the same linearization technique for F2" +1 as was done for Fl_, H 
the following system of equations in 64: 


M 


52 I<S = FSj 


1 = 1 

The variables K S and FS are defined as 
m m 


•£ks,, = Y. Ki - ~ eAx 


i=\ 


M 

-2 cr +1 E + cr +1 e K5 a 

i = l 


M 

L 

«=l 


(3.65) 

(3.66) 

(3.67) 

(3.68) 

(3.69) 
produces 


(3.70) 
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and 


M M 


M M 

+ cr' e e Kink k?m: + cr' e E Ki0 >* «?“: 

i=i jt=i * =1 fc=1 

M M 

- ^'EE^W 

1=1 Ar=l 


{ M 

-2[«cr +l +(i-«)cr]E A4 „ m. 

M 

[^cr +1 + (i - *)cr] xi ^5j. Mi 


+ 


i= l 
M M 


+ [0C 2 n+1 + (1 - 0)C2 n ] XI H # 8,-* fc6 * 63 ^ 

1 . = 1 fc= 1 

M Af 

[0C 3 n+1 + (1 - 0)C 3 "1 XI H /a0 i«* 66 * 63fc 

L i i 


+ 


;=i fc=i 
M M 


+ 


[0C 3 n+1 + (1 - o)cf\ XI XI K 9 i* b6ib3k 

1 i = 1 *=1 

M 

[«c; + ' +(i-«)c,"]E A6 J . il6 ' 


1 = 1 
Af 


Af "j 

+ [C 4 ” +1 + (i-«)cj]E/n„^|. 


(3.71) 


(3.72) 


Boundary conditions for r and 5 at the boundary-layer edge are automatically imposed 
through the use of the (1 - u) term in the trial solutions. As note previously, no wall 
boundary condition for r, and hence A61” +1 , is imposed. The boundary condition at the 
wall for s, and hence AM" +1 , is given by 

AM" +1 = Abl^s w . (3-73) 


The system of equations given by Equations 3.64 and 3.70 can be solved without 
iteration by the Thomas algorithm at each streamwise location using a variable step size. 
The variable step size is based on a comparison of the linearized and nonlinearized values 
of F 1" +1 and F2" +1 . The maximum relative error expression is given by 

Fl n+1 (linearized) - F 1 " +1 (nonlinearized) 

A — max Fl” +1 (nonlinearized) 

F2” +1 (linearized) - F2? +l (nonlinearized) ^ 3 74) 

F2" +1 (nonlinearized) 
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Alternately, the linearization error 
allowable ratio of 


A61 


n + 1 


61" 


and 


A64 


n + 1 


64? 


can be efficiently controlled by enforcing a maximum 
. The maximum ratio expression is given by 


A = max 


A61" +1 A64" +1 
61" ’ 64" 


(3.75) 


Given a maximum and minimum tolerance A m ,„ and A max , respectively, the step size is 
modified according to the following procedure: 


A mt „ < A < A mox step size is unchanged 
A > A mQX step size is halved 

A < A min step size is doubled 


3.4 Self-Similar Solutions 

The self-similar solution may be used near the leading edge to obtain the needed initial 
conditions. In order to solve Equation 3.64 for the shear stress at the (n + l) th streamwise 
location, it is necessary to know the shear stress at the n th streamwise location. A number 
of techniques have been used to obtain self-similar solutions for specific geometries within 
certain regions. The classic Falkner-Skan and Illingworth-Stewartson series can produce 
self-similar solutions for stagnation point and wedge flows with heat transfer [20], [21]. This 
same method may be applied to the traditional Dorodnitsyn equations of motion (laminar 
form of Equations 3.17, 3.20, and 3.22) to provide the initial r and 5 values. Even though 
the governing equations of interest are turbulent, the laminar self-similar solutions are valid 
in the region near the stagnation point where transition to turbulence has not yet occurred. 

At a small distance from the stagnation point, the velocity u e can be approximated by 

U e = UooC\ ( 3 - 76 ) 


where m is related to the pressure gradient coefficient /? as follows: 

2m 


/? = 


m + 1 


(3.77) 


Also, near the stagnation point A£ = £ - 0 and Ax = x - 0 such that Equation 3.13 can be 
written as 


U co Poo 


(3.78) 
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The similarity parameter a, 
by 


used to reduce the number of independent variables, is given 



(3.79) 


3.4.1 Momentum Similarity 

The development of the momentum similarity equation takes the same form as the classical 
development. Using the no-slip boundary condition and integrating with respect to rj, the 
continuity equation becomes 

w = f udrj. (3.80) 

ai Jo 

Representing the nondimensional velocity u as 

• ** =v, (3.8i) 


u = 


da 


and substituting into Equation 3.80 gives 


d 2 a 


d(.di) uj 0Lu/' 

w ~ ~ 72 ^ da ^ ' 


da 

di 


(3.82) 

(nr «" 

Substituting Equations 3.81 and 3.82 into the laminar form of Equation 3.20 generates a 
momentum similarity equation of the following form: 


1 du h 

\Il'" _L H ( 1 _ W' 2 ) _ didrtXY^" -- 0 

' • n ( da\ 


d 2 a 






(3.83) 


The coefficients are evaluated by obtaining the same form of similarity equation as White 
[20], which requires 




1 du e 
u € di 


and 


m 


Cofe) 3 , 


d 2 a 


Therefore, the final momentum similarity equation can be written as 

4-"' + 0 (1 - 'V' 2 ) + (^) VP" = 0, 

with the following boundary conditions: 

vp(0) = 'E'(O) = 0 and ^'(oo) = 1 


(3.84) 


(3.85) 


(3.86) 


(3.87) 
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3.4.2 Energy Similarity 

Representing the nondimensional enthalpy parameter s as 

5 = Cl(a), (3-88) 

and substituting into the laminar form of Equation 3.22 along with Equations 3.81 and 3.82 
generates an energy similarity equation of the following form. 

n" + - IPr/W'd + 2 ^ (1 - Pr) + W"] = 0 (3.89) 

m + 1 (7 + 1 ) 

The energy similarity coefficients are evaluated according to the previous definitions for 0 
and m. Heat-transfer boundary conditions are given by 

Q( 0 ) = = (l - \ and Q(oo) = 0, (3.90) 

V 1 0 / ti e 

where T w and T 0 denote wall and stagnation temperatures, respectively. If an adiabatic wall 
condition is prescribed, the boundary conditions are given as 

Q'( 0 ) = 0 and Q(oo) = 0. (3.91) 

Equations 3.86 and 3.89 define a two-point boundary-value problem (BVP). Using 
bisection for the shooting method along with a 4 th-order Runge-Kutta scheme, the BVP 
can be solved to produce an initial set of r and s values that can be used to begin the 
marching routine. The value of r is proportional to 4'" and is calculated by the following 

equation: , 

_ du _ da du _ ( (m + 1) Uoopoo V'P" (3.92) 

dr/ di i da \ 2Co u e p e x J 


3.5 Turbulence Models 


A turbulent boundary layer can be regarded approximately as a composite layer made up 
of inner and outer regions which arise due to the different response of the fluid to shear and 
pressure gradients near the wall [9]. In order to determine the eddy-viscosity p t in the inner 
and outer regions, two models based on the mixing-length concept are used in the SDG 
formulation. 

Prandtl proposed that each turbulent fluctuation could be related to a length scale and 


velocity gradient: 




pi 2 


du 

dy 


(3.93) 
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The mixing length l is related to the flow conditions and takes on different values in the 
inner and outer turbulent flow regions. 

A single inner region eddy-viscosity model, which follows the form used by Cebeci- 
Smith, was considered in this analysis [9]. The van Driest modified mixing-length formula 
in terms of the Dorodnitsyn transformation is given by 

Poo U e p^l 2 r, (3.94) 

H Co p e Ulc 

where / is the mixing length (k = 0.41) 

l = ny — e y+ /' 4+ ) . (3.95) 

The parameter A + is dependent on the pressure gradient according to 

4 + = r i ( 3 - 96 > 

(1 + 10 p+)s 


where A 0 = 26, 


Co Pe U(X) dp e 
\[Rt Poo p\, U \ dx ’ 


and 


( Co Pe r w \ * 

\ Poo Pw ) 


(3.97) 


The mixing length parameter / is related to the coordinates y + and y which are obtained 
from the solution by 


y + = 


Re Poo pi} 1 
Cq Pe Uc 


-y, where y = — f ==— / — (3.98) 
y ’ y \fRe u e Jo pr 


Two outer region eddy-viscosity formulations were considered in this analysis. The 
first outer region eddy- viscosity model considered is based on the form used by Cebeci and 
Smith [9]. The outer model is modified to account for intermittency near the edge of the 
boundary layer 6: 

fit 0.0168 Re poo u e ^ 399 ^ 

p, [1 + 5.5(y/<5) 6 ] Co p e u 0 o 

The second outer region eddy- viscosity model considered is based on the form developed 
by Baldwin and Lomax [22], Instead of using the boundary layer and displacement thickness 
as parameters in the outer formulation, the Baldwin-Lomax model uses certain maximum 
functions occurring in the boundary layer. The eddy-viscosity relation is given by 


[It 0.0168 Cep Re Poo ^max 

— = 7 — : “77 "77 P 2/max> 

fi [1 + 5.5(CklebJ//ymax) ] Cq Pe u oo 


(3.100) 
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where t/ max is the value of y corresponding to F max . The parameters Cep and Ckieb are 
calculated by the following expressions in order to fit the known properties of Coles wake 
law and equilibrium pressure gradients [23]: 


, 3 — 4 Ckieb 

,CP = 2C k |eb (2 - 3 C k l e b + C k 3 leb ) 


(3.101) 


2 0.01312 

Ckieb - 3 - 0.1724 + (i T 


2/max dUg 

u T dx 


(3.102) 


The choice between the inner and outer eddy-viscosity formula is made by taking the 
smaller value; the changeover typically occurs at u ~ 0.7. 
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Chapter 4 

Viscous-Inviscid Interaction 


In many flows of practical aerodynamic interest, the effects of viscosity and turbulence are 
confined to a relatively thin shear layer near the airfoil surface and wake. The flow over an 
airfoil can therefore be divided into two regions: an inner region, where viscous effects are 
important, and an outer inviscid region. Modelling the viscous and inviscid flow regions 
separately while providing a mechanism by which each solution influences the other is 
called Viscous-Inviscid Interaction (VII). A VII method organizes the viscous and inviscid 
parts of the overall solution to interact in an iterative way so that convergence of the final 
solution is achieved as economically and accurately as possible. The principle interaction 
between the regions arises from the displacement effect of the shear layers which leads to a 
thickening of the equivalent body with a corresponding change in the surface pressure [24]. 
The resulting interaction is labeled as either weak or strong, depending on the change in 
pressure and the extent to which higher-order viscous effects influence the overall solution. 
VII, based on a direct relationship between the viscous and inviscid regions of flow, is 
applicable as long as the disturbances to the inviscid flow due to the viscous displacement 
effect are small [25]. 

4.1 Viscous-Inviscid Coupling 

The classical interaction approach is to obtain an approximation to the inviscid flow, extract 
velocity and pressure from the inviscid solution, use the external conditions to obtain 
an approximation to the viscous flow, extract displacement thickness from the viscous 
solution, and use the viscous parameters to modify the original body geometry and obtain 
another estimate of the inviscid flow. An alternative to adding the displacement thickness 
distribution to the original body thickness is to impose a transpiration velocity boundary 
condition at the body surface. These iterative cycles, depicted in Figure 4.1, continue until 
the inviscid and viscous solutions are converged and compatible. 



31 



Figure 4.1 Viscous-Inviscid Interaction Scheme 


4.1.1 Displacement Thickness 

The displacement thickness is the height by which a streamline is displaced upward by the 
presence of the boundary layer. Considering the flow over a flat surface as depicted in 
Figure 4.2, the no-slip condition at the wall causes a partial obstruction of the freestream 
flow. This results in an upward deflection of the streamline passing through point h at 
station 1 by a distance 6 * at station 2. Because the flat surface and the streamline form 
the boundaries of a streamtube, the mass flow across stations 1 and 2 must be equal and is 

expressed by K K 

J p e u e dy = J pudy + p e u e fi m - (^. 1 ) 

Rearranging gives a familiar form of the compressible displacement thickness. 
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The nondimensionalized compressible displacement thickness in terms of the Dorodnitsyn 
formulation is given by 
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Figure 4.2 Viscous Flow Momentum Defect 


The advantage of using the displacement thickness as the coupling mechanism between 
the viscous and inviscid solutions is its simplicity. The disadvantage of the displacement 
thickness approach is that the inviscid grid must be generated after each viscous iteration 
because the body about which the inviscid flow is computed changes when the displacement 
thickness alters the effective body. An additional disadvantage of the displacement thickness 
approach is the possibility of supercritical VII which does not allow smooth transition to 
separation [26]. It has been shown that the classical displacement thickness interaction 
becomes supercritical when « 1.5 -» 2.0, depending on the history of the turbulent 
boundary layer [27]. 

4.1.2 Transpiration Velocity 

The transpiration velocity is an inviscid normal -velocity boundary condition which is 
imposed at the body surface to simulate the displacement of the inviscid flow by the viscous 
flow momentum defect [24], An expression for the transpiration velocity can be obtained 
by integrating the difference between the inviscid and viscous continuity equations across 
the boundary layer while applying the Prandtl boundary-layer and uniform inviscid-flow 
assumptions [26]. The transpiration velocity expression is given by 

v t = (peiie6 m ) ■ 

p e OX V ' 


(4.4) 



33 


The nondimensionalized compressible transpiration velocity in terms of the Dorodnitsyn 
formulation is given by 


d6‘ 

v t = u e — + 

dx 


( ft 

\ dx Pe dx j 


(4.5) 


The transpiration velocity distribution should result in a inviscid streamline coincident with 
the effective body obtained with the displacement thickness approach. The advantages of 
the transpiration velocity approach are that the inviscid grid need not be regenerated after 
each viscous iteration and that the interaction always allows smooth transition to separation 
[27]. 


4.2 Euler Equation Solver 

A fast and robust two-dimensional Euler code (GAUSS2), developed by Dr. Peter M. 
Hartwich of Vigyan Inc., is used as the inviscid flow solver in the VII scheme [28], [29]. 
The method uses a floating shock-fitting technique that has been combined with a second- 
order accurate upwind scheme based on the split-coefficient-matrix (SCM) method and 
with a time-implicit, diagonalized approximate-factorization (AF) algorithm. The result is 
a fast and robust two-dimensional Euler code that produces accurate solutions for shocked 
flows on crude meshes which are not adapted to the shock fronts. 

The equations for two-dimensional, compressible, and nonconservative Euler equations 
for a polytropic gas at constant 7 are given in general coordinates as 

Qv + AQ^ 4- BQv = 0, (4-6) 


where 

Q = (a',u',t/,s') T , 


(4.7) 


and the vectors A and B are defined by a coefficient matrix C [28]. The primes denote 
inviscid dependent variables where a', u\ v\ and s' are speed of sound, cartesian velocities, 
and entropy, respectively. The dependent variables have been nondimensionalized as 
follows: 


x 


t 





_P_ 

Poo 


(4.8) 
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4.2.1 Viscous-Inviscid Interfaces 

In order to allow the boundary-layer and Euler methods to work together in an VII scheme, 
it is necessary to develop a set of expressions which convert the variables from one type 
of nondimensionalization to another. The boundary-layer method employs a body-normal 
coordinate system while the Euler method employs a cartesian coordinate system, so a 
coordinate system conversion must be performed. Also, the methods use different reference 
values to nondimensionalize velocity and pressure, so a dependent variable conversion must 
be performed. The variable conversions which account for both the coordinate system and 
nondimensionalization differences are summarized in Table 4.1, where 6 is the local body 
angle. 


Table 4.1 Viscous-Inviscid Variable Conversion 


Variable 

Euler — ► Boundary Layer 

Length 

x = \Jx a + y a 

Velocity 

u = (u'cos0 -1- u'sin 9)y/pZ 

Pressure 

11 

8 
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Chapter 5 


Numerical Results 


5.1 Convergence Properties 


The convergence properties of the Dorodnitsyn finite element approximation of r are plotted 
in Figure 5.1. The average discrete root-mean-square of the relative error is calculated by 
the following equation: 


Trms 




(M-\) 


(5.1) 


The exact value of r is obtained from a solution of the similarity equations developed in 
Section 3.4.1 for flow over a flat plate with zero pressure gradient (/? = 0). The 4th-order 
Runge-Kutta step size was set to Aa = 0.01. The results presented in Figure 5.1 are for 
8 values of Au corresponding to a discretization of 3, 4, 5, 7, 9, 11, 15, and 19 nodes at 
x = 0.5 along the plate. In order to evaluate the convergence properties of the spatial 
discretization across the boundary layer, the step size was held at a small constant value to 
minimize the influence of the marching routine. Examination of the figure indicates that 
the use of linear elements achieves the theoretically expected second-order convergence. 
It can be seen that the convergence rate of the scheme decreases slightly as the number 
of nodes increases. The decreasing convergence rate can be attributed to the decrease 
in accuracy of the linearized marching scheme with increasing number of nodes. An 
extensive analysis of convergence properties for favorable and adverse pressure-gradient 
cases as well as quadratic elements may be found in Reference [18]. Fletcher and Fleet 
showed that the accuracy of linear elements on coarse grids are comparable to quadratic 
elements, and the convergence properties for the adverse and favorable pressure-gradient 
cases are comparable to the zero- gradient case. 


5.2 Flat Plate 

An assessment of the accuracy of the semi-discrete Galerkin (SDG) method may be per- 
formed by examining the ability of the method to reproduce the well-known compressible 
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Figure 5.1 Nondimensional Shear Stress Convergence 
Results for a Flat Plate at Moo = 0.800, Re = 5000 


laminar boundary-layer solution over a thermally-insulated flat plate with zero pressure 
gradient. The computations were performed at a Mach number of 0.8 and a Reynolds 
number based on plate length of 5000. A discretization of 9 nodes was used across 
the boundary layer. The marching routine step size was efficiently controlled by setting 
A m , n = 1.0 x 10" 4 = 0.lA mar to maintain the ratio of -rf*- and -^p. Again, an exact 
analytical solution for this flow may be obtained by solution of the Falkner-Skan similarity 
equation [20]. Figure 5.2 shows excellent agreement between the exact and computed 
nondimensional shear-stress values at x = 0.5. The similarity property of the solution was 
also verified by examining the profiles at different stations along the length of the plate. 
Excellent agreement was observed except for stations close to the leading edge where effects 
of the stagnation-point flow are present. The skin friction along the plate was calculated 

using the relation 

C, = T. . (5.2) 

V Re ul, poo 

To illustrate the high accuracy of the Dorodnitsyn formulation on extremely rough grids, 
both 3 and 9-node coefficient of friction solutions are plotted along with the exact solution 
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Figure 5.2 Comparison of Computed and Exact Nondimensional Shear Stress 
Profiles for a Flat Plate at x = 0.5, Moo — 0.800, Re = 5000 


in Figure 5.3. There is excellent agreement between the 9-node and exact solutions. The 
3-node solution produced good results except in the region very close to the leading edge. 

5.3 NACA 0012 Airfoil 

The NACA 0012 was chosen as the primary test airfoil for the compressible turbulent VII 
validation. The NACA 0012 is a symmetric airfoil which has been tested both computa- 
tionally and experimentally by a great number of researchers. Of particular importance to 
the VII scheme is that the experimental data at zero angle of attack is not affected by wall 
interference due to lift. The airfoil is tested for attached turbulent flow at high Reynolds 
numbers and transonic speeds as shown in Table 5.1. In order to minimize the influence 
of wind tunnel wall-interference effects in the comparison of experimental and computed 
results, the numerical angle of attack was varied in each case to match the computed lift 
and the experimental normal force coefficients. In each of the NACA 0012 cases, the 
161x31 C-grid shown in Figure 5.4 was used by GAUSS2 to solve the inviscid equations of 
motion. The outer boundary of the C-grid is located 5 to 6 chord lengths from the airfoil in 
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Figure 5.3 3 and 9-Node Comparison of Computed and Exact Skin 

Friction for a Flat Plate at M oo = 0.800, Re = 5000 


all directions. The transpiration velocity boundary condition was used as the VII coupling 
mechanism and was relaxed after each global VII iteration according to 

vl el = V q t +(\-0 J )vr\ ( 5 - 3 ) 

where q denotes a particular global VII iteration and u is a relaxation parameter. The 
relaxation parameter was set equal to 0.1 for all NACA 0012 calculations. Convergence 
of the scheme was assumed when lift and total drag coefficients changed less than 0.1% 
between global iterations. The maximum number of global VII iterations was set at 50, 
and the ratio of local inviscid to viscous iterations was set at 100. A boundary-layer 
discretization of 9 nodes was used and the marching routine was controlled in the same 
manner as the flat-plate calculations by setting A = 5.0 x 10 -4 . 

53.1 NACA 0012 - Case A 

This test case consists of a NACA 0012 airfoil at a numerical angle of attack of -0.14, a 
freestream Mach number of 0.499, and a Reynolds number of 9.0 x 10 6 . The flow is attached 
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Table 5.1 NACA 0012 Test Cases 


Case 

Moo 

Re 

&exp 

&num 

Ref. 

A 

0.499 

9.0 x 10 6 

-0.14 

-0.14 

[30] 

B 

0.700 

9.0 x 10 6 

1.86 

1.37 

[30] 


and subsonic over the entire airfoil surface. Transition to turbulence was numerically tripped 
at the leading edge. The computed coefficient of pressure compares reasonably well with 
experiment but is slightly underpredicted over the aft section to the trailing edge of the 
airfoil as shown in Figure 5.5. No experimental data was readily available to compare the 
computed friction coefficient and displacement thickness, but both Cj and <5 are plotted in 
Figures 5.6 and 5.7 for future reference. 

53.2 NACA 0012 - Case B 

This test case consists of a NACA 0012 airfoil at a numerical angle of attack of 1.37, a 
freestream Mach number of 0.700, and a Reynolds number of 9.0 x 10 6 . For this case, 
the flow is attached and just slightly supersonic near the leading-edge upper surface. A 
comparison of experimental and calculated coefficient of pressure is shown in Figure 5.8. 
The predicted C v compares well with experiment but slightly overpredicts the peak value 
at the supersonic region near the leading edge and deviates at the trailing edge. Again, the 
computed skin friction and displacement thickness values are plotted in Figures 5.9 and 
5.10 for future reference. 

5.3 3 Aerodynamic Characteristics 

Figure 5.11 shows a comparison of lift coefficient vs. angle of attack for the NACA 
0012 airfoil at = 0 700 and Re = 9.0 x 10 6 . The SDG VII scheme overpredicts 
coefficient of lift at low angles of attack, but the prediction improves at higher angles of 
attack. For angles of attack above 1.54, the flow is more strongly transonic and eventually 
separates. Drag polar comparisons are displayed in Figure 5.12 for the NACA 0012 airfoil 
at Moo = 0.700 and Re = 9.0 x 10 6 . The data again shows that the computed lift coefficient 
is overpredicted at lower angles of attack. It is also noted that the computed drag coefficient 
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Figure 5.5 Comparison of Pressure Coefficient Distribution for the NACA 
0012 Airfoil at Moo = 0.499, a num = -0.14, Re = 9.0 x 10 6 



Figure 5.6 Skin Friction Coefficient Distribution for the NACA 0012 Airfoil 
Upper-Surface at = 0.499, a num = -0.14, Re = 9.0 x 10 6 
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Figure 5.9 Skin Friction Coefficient Distribution for the NACA 0012 Airfoil 
Upper-Surface at Moo = 0.700, a num = 1.37, Re = 9.0 x 10 6 



Figure 5.10 Upper-Surface Displacement Thickness for the NACA 
0012 Airfoil at = 0.700, a num = 1.37, Re = 9.0 x 10 6 
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is underpredicted at lower angles of attack with an improved prediction at higher incidence. 
The flowfield is subsonic for C L values of approximately 0.2 and lower. Drag values below 
this point correspond to pressure plus skin friction drag. Drag values above this point have, 
in addition, a wave drag component. Transonic drag-rise characteristics for the NACA 
0012 airfoil at zero-lift conditions are displayed in Figure 5.13. The turbulent boundary 
layer was numerically tripped at the leading edge, and all computations were performed at 
a Reynolds number of 9 million. The cross-hatch range of experimental values is based 
on a “best of six” set of data as described by McCroskey [32]. An underprediction of the 
drag coefficient is observed at lower Mach numbers with a trend toward improved results 
at higher Mach numbers. 


0.8 


0.6 


Cl 0.4 


0.2 


0 


Figure 5.11 Comparison of Lift Coefficient vs. Angle of Attack for 
the NACA 0012 Airfoil at Moo = 0.700, Re = 9.0 x 10 6 


5.4 RAE 2822 Airfoil 

The RAE 2822 was chosen as a supplementary test airfoil for the compressible turbulent 
VII interaction. The RAE 2822 is a supercritical airfoil with a moderate amount of aft 
camber which poses a challenge in achieving VII convergence. Also, the experimental 
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Figure 5.12 Comparison of Lift vs. Drag Polars for the 
NACA 0012 Airfoil at M = 0.700, Re = 9.0 x 10 6 



Figure 5.13 Comparison of Computed and Measured Transonic Drag-Rise 
Characteristics for the NACA 0012 Airfoil at o B um = 0* = 9.0 x 10 
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data in Reference [31] contains a number of boundary-layer and wake parameters such 
as displacement thickness, momentum thickness, and skin friction which are helpful in 
validating the boundary-layer solution. The airfoil is tested for attached turbulent flow at 
high Reynolds numbers and transonic speeds as shown in Table 5.2. A C-grid similar to the 


Table 5.2 RAE 2822 Test Cases 


Case 

Moo 

Re 

&exp 

^num 

Ref. 

A 

0.676 

5.7 x 10 6 

-2.18 

-2.18 

[31] 

B 

0.676 

5.7 x 10 6 

2.40 

1.90 

[31] 

C 

0.725 

6.5 x 10 6 

2.55 

2.10 

[31] 


one use for the NACA 0012 test cases was used for all of the RAE 2822 calculations. The 
relaxation parameter was reduced to u> = 0.05 which produced a slower but more stable VII 
convergence. Consequently, the maximum number of global VII iterations was increased 
to 75. The discretization and marching-routine controls were set at values equal to those 
used in the NACA 0012 test cases. For the all of the RAE 2882 cases, the skin friction 
coefficient is based on the boundary-layer edge dynamic pressure. 

5.4.1 RAE 2822 - Case A 

This test case consists of a RAE 2822 airfoil at a numerical angle of attack of -2.18, a 
freestream Mach number of 0.676, and a Reynolds number of 5.7 x 10 . The flow in 
this case is subsonic and attached over the entire airfoil surface. The computed coefficient 
of pressure shown in Figure 5.14 compares well with the experimental values. The C p 
is slightly underpredicted on the lower surface near the leading edge and again deviates 
from experiment at the trailing edge. The coefficient of friction and displacement thickness 
results given in Figures 5.15 and 5.16 have good agreement with the experimental values. 
As seen in Figure 5.15, turbulence was tripped at x = 0.11. Figure 5.17 shows the 
velocity profiles in the boundary layer at three x-locations. The velocity profile deviation 
from experimental values is attributed to the fact that the SDG method calculates velocity 
indirectly since shear stress is a dependent variable. 
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Figure 5.14 Comparison of Pressure Coefficient Distribution for the RAE 
2822 Airfoil at M<* = 0.676, o num = -2.18, Re = 5.7 x 10 6 



Figure 5.15 Comparison of Upper-Surface Skin Friction Coefficient for the RAE 
2822 Airfoil at M <*, = 0.676, a num = —2.18, Re = 5.7 x 10 6 
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Figure 5.16 Upper-Surface Displacement Thickness Comparison for the RAE 
2822 Airfoil at Moo = 0.676, a num = -2.18, Re = 5.7 x 10 6 



Figure 5.17 Comparison of Computed and Measured Velocity Profiles for the 
RAE 2822 Airfoil at = 0.676, a num = -2.18, Re = 5.7 x 10 6 
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5.4.2 RAE 2822 - Case B 

This test case consists of a RAE 2822 airfoil at a numerical angle of attack of 1.90, a 
freestream Mach number of 0.676, and a Reynolds number of 5.7 x 10 6 . The flow is 
attached over the entire surface and slightly supersonic near the leading edge on the upper 
surface of the airfoil. Transition to turbulence was numerically tripped at i = 0.03. The 
computed coefficient of pressure is compared with experimental values in Figure 5.18. The 
C p prediction on the lower surface of the airfoil is in good agreement with experiment except 
at the trailing edge, and the upper-surface prediction is slightly below experimental values 
over the aft section of the airfoil. It should be noted that the pressure coefficient prediction 
in the supersonic flow region is in good agreement with experiment. Both coefficient of 
friction and displacement thickness values compare well with the experimental values as 

shown in Figures 5.19-5.20. 



Figure 5.18 Comparison of Pressure Coefficient Distribution for the RAE 
2822 Airfoil at Moo = 0.676, a num = 1.90, Re = 5.7 x 10 6 




SDG VII 
Experiment 



Figure 5.19 Comparison of Upper-Surface Skin Friction Coefficient for the 
RAE 2822 Airfoil at Moo = 0.676, a num = 1-90, Re = 5.7 x 10 6 



Figure 5.20 Upper-Surface Displacement Thickness Comparison for the RAE 
2822 Airfoil at = 0.676, a num = 1.90, Re = 5.7 x 10 6 
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5.4.3 RAE 2822 - Case C 

This test case consists of a RAE 2822 airfoil at a numerical angle of attack of 2.10, a 
ffeestream Mach number of 0.725, and a Reynolds number of 6.5 x 10 6 . The transition to 
turbulence was again numerically tripped at x = 0.03. For this case, the flow is attached 
and supersonic on the upper surface where a moderately strong shock wave is experienced. 
A comparison of experimental and calculated coefficient of pressure is shown in Figure 
5.21. The predicted C p compares relatively well with experiment except at the shock wave 
which is predicted upstream of the experimental result. The coefficient of friction and 
displacement thickness results given in Figures 5.22 and 5.23 have good agreement with 
the experimental values. 



Figure 5.21 Comparison of Pressure Coefficient Distribution for the RAE 
2822 Airfoil at Moo = 0.725, a num = 2.10, Re = 6.5 x 10 6 


5.5 Grid Refinement Study 

A grid refinement study was performed to demonstrate the sensitivity of computed force 
coefficients to the grid spacing used in the SDG VII scheme. The NACA 0012 airfoil 
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Figure 5.22 Comparison of Upper-Surface Skin Friction Coefficient for the 
RAE 2822 Airfoil at = 0.725, a num = 2.10, Re = 6.5 x 10 6 



X 


Figure 5.23 Upper-Surface Displacement Thickness Comparison for the RAE 
2822 Airfoil at = 0.725, a num = 2.10, Re = 6.5 x 10 6 
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was solved at flow conditions given in Case A for a number of viscous finite element 
discretizations and inviscid grid sizes. Figure 5.24 is a plot of computed lift coefficient 
vs. finite element space width for coarse (81 x 17) and fine (161 x 33) inviscid grids. 
The trend is to decrease the lift coefficient with increasing space width. The effect of grid 
spacing on the computed drag coefficient is shown in Figure 5.25, where the trend is to 
increase the drag coefficient with increasing space width. The relatively small slope of 
each curve indicates the method produces reasonable drag levels on coarse grids, which is 
a highly-desirable characteristic of any computational method. 



Figure 5.24 Comparison of Lift Coefficient vs. Boundary-Layer 
Space Width for the N ACA 00 1 2 Airfoil - Case A 


5.6 Aerodynamic Force Coefficients 

A summary of computed lift and drag coefficients for the SDG VII method is displayed 
in Table 5.3. All drag coefficient values are in terms of drag counts where C D = 1 count 
is equivalent to Cd = 0.0001. Since the numerical angle of attack was varied in each 
case to match computed lift and experimental normal force coefficients, it is important to 
note that a num was within the accepted corrected angle of attack range for wind tunnel 
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Figure 5.25 Comparison of Drag Coefficient vs. Average 
Grid Spacing for the NACA 0012 Airfoil - Case A 


wall-interference effects [33]. The computed coefficient of drag results for the SDG VII 
scheme were on average predicted approximately 15% lower than experimentally obtained 
results. The underprediction of drag is more pronounced at lower angle of attack cases. 
The current state-of-the-art capabilities for attached flow lift and drag predictions are ±3% 
range for lift and ±5% range for drag [33]. 




Table 5.3 Force Coefficient Comparison 


NACA 0012 

Numerical 

Experimental 

Case 

C L 

Cd 

Cpp 

Cdv 

Cl 

C D 

A 

-0.019 

58 

-7 

65 

-0.013 

11 

B 

0.242 

70 

7 

63 

0.241 

79 

RAE 2822 

Numerical 

Experimental 

Case 

Ci ^ 

Cp n 

Cop 

Cdv 

Cl 

C D 

A 

-0.121 

70 

7 

63 

-0.115 

79 

B 

0.565 

75 

8 

67 

0.566 

85 

C 

0.658 

86 

21 

65 

0.658 

107 
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Chapter 6 
Conclusions 


The application of the Dorodnitysn transformations to the boundary-layer equations pro- 
vides a spatial coordinate which automatically follows boundary-layer growth and gives 
high resolution near the wall which is important in turbulent flows where near- wall velocity 
gradients are large. The transformed spatial domain is effectively modelled with finite ele- 
ments to provide accurate results on coarse grids. Modelling the group terms for viscosity 
by the group finite element method leads to an important economy in the formulation since 
fit need only be evaluated at the nodes. The use of a non-iterative streamwise marching 
routine also adds computational economy to the method. The resulting boundary layer 
code can be used to solve any smoothly-connected airfoil shape without modification and 
is easily coupled with existing inviscid flow solvers. The transpiration velocity approach 
used to couple the particularly fast and accurate Euler solver used in this analysis serves as 
a computationally efficient and easily implemented VII coupling mechanism. Convergence 
of the overall interaction procedure for both the NACA 0012 and RAE 2822 airfoils was 
achieved in relatively few global iterations while achieving results of adequate engineering 
accuracy. 

Future work with the SDG method should include the use of a maximum reversed 
flow velocity concept in order to successfully model separation while still retaining a finite 
spatial grid. Also, the marching routine could be further improved by possibly using a 
trigonometric streamwise coordinate transformation to smooth external velocity gradients 
at the leading and trailing edges of the airfoil. 
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Appendix A 


Galerkin Integrals 
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Introduction 

A method is developed for solution of the two-dimensional, steady, compressible, turbulent 
boundary-layer equations and is coupled to an existing Euler solver for attached transonic 
airfoil analysis problems. The boundary-layer formulation utilizes the semi-discrete Galerkin 
(SDG) method to model the spatial variables with linear finite elements and the time-like 
variables with finite differences. A Dorodnitsyn transformed system of equations is used to 
replace the infinite spatial domain with a finite domain thereby permitting the use of a uniform 
finite element grid which provides high resolution near the wall and automatically follows 
boundary-layer growth. The second-order accurate Crank-Nicholson scheme is applied along 
with a linearization method to take advantage of the parabolic nature of the boundary-layer 
equations and generate a non-iterative marching routine. The SDG code can be applied to any 
smoothly-connected airfoil shape without modification and can be coupled to any inviscid flow 
solver. A direct viscous-inviscid interaction is accomplished between the Euler (GAUSS2) and 
boundary-layer (SDGM) codes through the application of a transpiration velocity boundary 
condition. Gross effects of turbulence in the boundary layer are modelled through the use of a 
zero-equation algebraic Cebici-Smith or Baldwin-Lomax eddy-viscosity model. 

Program Structure 

There are three programs which comprise the current airfoil analysis package. The 
boundary-layer and Euler codes are called SDGM and GAUSS2, respectively. The interaction 
of the inner and outer region flow solvers is accomplished by a third control code called VII. 
There are a number of input and output date files which are utilized by the three codes. The 
control code VII is the only code which is explicitly invoked from the command line, and VII 
consequently uses internal system calls to execute either GAUSS2 or SDGM to iteratively 
perform the viscous and inviscid flow calculations. GAUSS2 utilizes a cartesian C-grid to 
solve for the entire inviscid flowfield while SDGM uses a body-normal coordinate system to 
solve for the viscous flowfield on the upper and lower surfaces of the airfoil. 


Data Files 

The following figure shows all of the data files which are either read or written by the three 
codes. Arrows indicate whether the data is being read or written by the codes and indicates the 
interaction of data files between the codes. Many of the data files are used for interim data 
storage purposes and are not of particular interest to the end-user of the programs, but they are 
included here for completeness. The following is a list of data files which should be properly 
edited and available to the codes in order to start a flow analysis procedure: 

Filename Description 

1) viinput.dat input data used by VII to control overall viscous inviscid interaction 

process 

2) sdgminputl .dat input data used by SDGM to control viscous flow calculations on 

upper surface of airfoil 



3) sdgminput2.dat input data used by SDGM to control viscous flow calculations on 

lower surface of airfoil 

4) siminputl.dat input data used by SDGM to control similarity solution calculations 

on upper surface of airfoil 

5) siminput2.dat input data used by SDGM to control similarity solution calculations 

on lower surface of airfoil 

6) input input data used by GAUSS2 to control inviscid flow calculations 

7) AIRGEO input data used by GAUSS2 which contains inviscid C-grid 

The following is a list of data files which contain the final output: 


Filename Description 


8) convergence.dat 

9) fcoeff.dat 

10) bl2ext.dat 

11) pcoeff.dat 


output data generated by VII which contains convergence history of 
aerodynamic force coefficients 

output data generated by SDGM which contains skin-friction 
coefficient over the airfoil surface at streamwise locations based on 
the viscous marching algorithm 

output data generated by SDGM which contains extrapolated 
integral quantities used in interaction algorithms at streamwise 
locations based on the inviscid grid 

output data generated by GAUSS2 which contains pressure 
coefficient over the airfoil surface at streamwise locations based on 
the inviscid grid 


Structure and Guidelines 

The following is further definition of the necessary data files with an example format, variable 
definitions, and parameter guidelines: 

1) viinput.dat 

iinit 

istart 

iend 

abschangel 

abschange2 

iinit - arbitrary starting number for global vii iterations ( typically 1 ) . 

istart - number of iterations to wait before checking for convergence ( typically 1/relax ) 
iend - maximum allowable number of global vii iterations ( typically 50 ) 
abschangel - minimum absolute change in lift coefficient to signal vii convergence 
( typically 0.001 ) 

abschange2 - minimum absolute change in total drag coefficient to signal vn convergence 
( typically 0.0001 ) 



2) sdgminput 1 .dat 

node 

theta 

P r 

re 

initx 

totx 

delx 

transx 

maxerror 

extrbl 

seprbl 

relax 

node - number of nodes in finite -element discretization ( typically betwen 7 and 13 ) 

theta - implicitness factor for crank-nicholson ( typically 0.5 ) 

pr - prandtl number ( typically 1.0) 

re - freestream reynolds number based on chord 

initx - starting cartesian x/c location ( typically 5.0 E-4 ) 

totx - ending cartesian x/c location ( typically 1.0 ) 

delx - initial finite-difference marching step ( typically 1.0 E-6 ) 

transx - percent of cartesian x/c that transition to turbulence occurs 

maxerror - maximum allowable ratio of the inverse shear stress nodal value and the 

change in the inverse shear stress nodal value ( typically 1 .0/(2 A (node+2)) ) 
extrbl - maximum value of the inverse shear stress nodal value that indicates the onset of 
separation and the need to extrapolate future integral values ( typically 10.0 ) 
seprbl - maximum value of the inverse shear stress nodal value that indicates separation 
( typically 100.0 ) 

relax - vii relaxation parameter for the coupling mechanism ( typically 0.1 ) 

3) sdgminput2.dat - same as above 

4) siminputl.dat 

alphastart alphaend h eps 

zguess(l) zguess(2) zguess(3) zguess(4) 

ebctype 

twto 

m 

alphastart - starting value of similarity parameter alpha ( typically 0 ) 

alphaend - asymptotic ending value of similarity parameter alpha ( typically 10 ) 

h - space width for similarity parameter alpha ( typically 0.01 ) 

eps - error tolerance for shooting method convergence ( typically 1.0 E-5 ) 

zguess(l) = lower bound for unknown momentum boundary condition ( typically 1.93 
for stagnation point flow ) 



zguess(2) = upper bound for unknown momentum boundary condition ( typically 1.94 
for stagnation point flow ) 

zguess(3) = lower bound for unknown energy boundary condition ( typically -10*twto 
for adiabatic case ) 

zguess(4) = upper bound for unknown energy boundary condition ( typically 10*twto 
for adiabatic case ) 

ebctype - flag to indicate type of energy wall-boundary condition ( 0=adiabatic, l=heat 
transfer ) 

twto - ratio of wall and stagnation temperatures ( typically 1.0 ) 
m - falkner-skan pressure gradient factor where beta = (2*m)/(m+l) ( typically 1.0 
for stagnation point airfoil flows) 

note: the most computationally efficient case is to set twto and pr = 1 

5) siminput2.dat - same as above 

6) input 
alpha mach 

cfl cflb ncyc iflaga iflagb sigma theta 
connect darcy xpl xp2 xp3 xp4 
hedge mose xmax 
ncycr 

alpha = angle of attack in degrees 
mach = freestream mach number 
cfl = ( typically 10.0 ) 

cflb = local courant number ( typically 1.0 E-3 ) 

ncyc = number of total cycles (typically > 500 ) 

iflaga = number of total cycles before shock fitting ( typically 50 ) 

iflagb = ( typically 9999 ) 

sigma = ( typically 1.015 ) 

theta = ( typically 60.0 ) 

connect = ( typically 0 ) 

darcy = ( typically 0.0 ) 

xpl = ( typically 1.0 ) 

xp2 = ( typically 1.0 ) 

xp3 = ( typically 1.0 ) 

xp4 = ( typically 1.0) .... . f , 

hedge = number of inviscid nodes from edge of inviscid grid to trailing edge of the 

airfoil ( typically 14 ) 
mose = ( typically 1.58 E-2 ) 
xmax = ( typically 7.0 ) 

ncycr = number of restart cycles ( typically < 100 ) 


7) 


AIRGEO - see example file 




